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ABSTRACT 

RANS  computations  around  airfoils  are  presented  for  flows  where  transition  takes  place  across  a  laminar 
separation  bubble.  Transition  locations  are  found  employing  the  eN -method  and  a  linear  stability  solver  that 
uses  the  velocity  and  temperature  profiles  from  the  RANS-solution.  Hereby ,  the  known  eN -method  is  applied 
for  the  hows  with  steady  onset  conditions ,  whereas  a  new  temporal  scheme  for  the  eN  -method  is  used  for 
the  flow  around  a  plunging  airfoil.  Obtained  numerical  results  are  compared  to  hotwire-,  pressure-  and  PIV- 
measurements  of  laminar  separation  bubbles. 

1  INTRODUCTION 

The  aerodynamics  of  low-speed  low-Reynolds-number  flows  over  airfoils  and  wings  has  been  studied  for 
decades  in  experimental  and  numerical  works.  However,  the  research  was  often  of  an  academic  nature  since 
there  were  few  technical  applications.  One  exception  was  the  design  of  model  airplanes  for  which  a  rich 
experimental  data  base  grew  over  the  years.  Here  the  pioneering  work  of  F.  Schmitz  [1]  is  mentioned,  and 
much  work  was  later  performed  by  researchers  at  Stuttgart  University,  Germany  [2],  [3].  More  recently,  M. 
Selig  and  co-workers  continued  to  develop  advanced  aerodynamic  design  methodologies  for  model  airplane 
airfoils  and  performed  a  large  amount  of  measurements  [4],  [  ].  The  subject  of  low-Reynolds  number  flows 
does  receive  more  attention  presently  as  advances  in  the  field  of  micro  system  technologies  enable  the  de¬ 
velopment  of  Micro  Aerial  Vehicles  (MAV)  with  a  broad  variety  of  applications.  In  this  context  it  is  useful 
to  view  the  range  of  existing  and  future  unmanned  aerial  vehicles  by  using  the  Reynolds  number  as  the  flow 
parameter,  see  Tab.  1.  Dimensional  analysis  states  that  the  Reynolds  number  is  the  governing  parameter  that 
determines  lift  and  drag  of  a  given  geometrical  shape  and  a  fixed  disturbance  level  of  the  incoming  freestream 
flow. 


MAV  long  term 

MAV  medium  term 

Current  MAV 

Small  drones 

Drones 

Mass 

10* 

100g 

5kg 

50  kg 

Reynolds  number 

103  - 104 

104  —  5  •  104 

5  •  104  —  2  - 105 

2  •  105  -  5  •  105 

5  •  105  —  2  - 106 

Wing  configuration 

Flapping 

Fixed,  rotary,  flapping 

Fixed,  rotary 

Fixed 

Fixed 

Tab.  1 :  Range  of  airfoil  Reynolds  numbers  relative  to  Micro  Air  Vehicle  weight 
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RANS  Simulation  of  the  Transitional  Flow  Around  Airfoils  at 
Low  Reynolds  Numbers  for  Steady  and  Unsteady  Onset  Conditions 

Very  small  aerial  vehicles  foreseen  in  the  long  term  will  most  likely  exhibit  a  laminar  flow.  At  larger 
Reynolds  numbers,  mixed  laminar  and  turbulent  flows  occur  with  transition  in  between.  Moreover,  it  turns  out 
that  aerodynamic  designs  for  Reynolds  numbers  in  the  range  between  104  and  106  face  a  flow  phenomenon 
called  laminar  separation  bubble  (LSB).  The  LSB  is  a  flow  where  a  laminar  separation  takes  place  that  is  in 
most  cases  caused  by  an  adverse  pressure  gradient  along  the  smooth  aerodynamic  surface.  Small  disturbances 
present  in  the  laminar  flow  are  strongly  amplified  in  the  shear  layer  of  the  separated  flow  and  rapid  transition 
to  turbulence  takes  place.  The  turbulence,  in  turn,  creates  a  large  momentum  transport  normal  to  the  shear 
layer  so  that  the  flow  reattaches  to  the  surface  and  a  closed  bubble  is  formed  in  the  time-averaged  mean. 

Laminar  separation  bubbles  can  have  large,  adverse  aerodynamic  effects.  Usually,  they  create  additional 
drag  as  they  displace  the  outer  inviscid  flow.  This  results  in  reduced  suction  over  the  forward  portion  of 
airfoils  and  wings  and  decreases  pressure  recovery  in  rear  parts,  so  the  additional  drag  is  mostly  pressure 
drag.  The  increase  of  pressure  drag  depends  on  the  size  of  the  LSB,  in  particular  on  its  thickness  in  the 
normal-wall  direction.  A  more  dramatic  effect  occurs  if  the  transition  process  in  the  separated  shear  layer 
is  relatively  slow  and  the  adverse  pressure  gradient  is  strong.  Then  turbulent  momentum  transport  is  not 
sufficient  to  close  the  bubble  and  a  large  separation  occurs  that  extends  right  to  the  trailing  edge.  This  causes 
a  sudden  loss  of  lift  and  a  strong  increase  of  drag  along  with  significant  hysteresis  effects  of  force  coefficients 
with  varying  angle  of  attack.  Note  that  the  break  up  of  LSBs  may  be  experienced  over  a  broad  range  of 
Reynolds  numbers. 


exponential  break  down 

receptivity  growth  nonlinear  to  turbulence 


laminar 
flow 


turbulent 

flow 


reversed  flow 


Fig.  1 :  Instability  and  transition  in  a  laminar  separation  bubble,  with  a  sketch  of  time-averaged  streamlines 


The  principal  flow  behaviour  of  a  laminar  separation  bubble  is  sketched  in  Fig.  1.  External  distortions, 
for  example  acoustic  waves  or  free  stream  turbulence,  generate  small,  harmonic  waves  within  the  laminar 
boundary  layer  upstream  of  the  separation.  The  behaviour  of  the  waves  of  the  laminar  boundary  layer  can  be 
described  by  linear  stability  theory  (LST)  within  a  certain  spatial  region:  unstable  waves  grow  exponentially 
while  travelling  downstream.  This  linear  process  constitutes  the  second  stage  of  transition.  In  the  third  stage 
of  transition  the  distortions  become  so  large  that  saturation  occurs  and  secondary  instabilities  can  grow  on 
the  distorted  boundary  layer  flow:  This  behaviour  is  denoted  with  nonlinear  interactions.  Finally,  the  number 
of  spatial  and  temporal  modes  grows  rapidly  and  the  ordered  laminar  structures  break  down  to  turbulence. 

Transition  in  2D  laminar  separation  bubbles  is  dominated  by  the  behaviour  of  the  primary  2D  disturbances 
within  the  laminar  boundary  layer  (called  Tollmien-Schlichting  modes),  according  to  the  research  of  Rist, 
Wtirz  and  Wagner  ([6],  [7]  ,  [8]).  Initial  3D-disturbances  are  also  amplified  but  numerical  and  experimental 
evidence  indicates  that  they  do  not  govern  the  transition  process  and  its  location. 

2  FLOW  MODELLING  FOR  ENGINEERING  COMPUTATIONS 

Aerodynamic  design  and  analysis  tools  used  for  flows  with  laminar  separation  bubbles  are  presently  restricted 
to  2D  steady  flows  [9]  and  they  are  rather  not  reliable  in  their  predictions  of  the  break-up  process  of  LSBs. 
For  the  future  design  of  MAV's  more  general  analysis  tools  are  sought,  which  can  be  applied  to  2D  airfoils 
as  well  as  3D  flows  over  more  complex  configurations.  Direct  Numerical  Simulation  (DNS)  of  the  flow 
during  aerodynamic  design  cycles  is  not  feasible  because  of  the  extreme  demands  on  computational  resources 
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associated  with  this  approach.  This  paper  therefore  presents  the  current  status  in  modelling  and  computation 
of  low-Reynolds  number  flows  with  laminar  separation  bubbles  using  the  Reynolds-averaged  Navier-Stokes 
equations.  This  engineering  approach  is  expected  to  complement  the  existing  approximate  methods  in  the 
near  future.  The  main  issues  using  this  approach  are  the  modelling  of  the  transition  laminar  turbulent  as  well 
as  accurate  turbulence  modelling. 


2.1  Transition  Modelling 

2.1.1  Linear  Stability  Theory 

For  2D  flows  the  Linear  Stability  Theory  uses  a  coordinate  system  (x,z)  that  is  aligned  with  the  local  edge 
velocity  of  the  undisturbed  laminar  mean  flow  such  that  x  is  streamwise  and  z  is  wall  normal.  The  baseflow  at 
a  given  x  position  is  represented  by  (U,0),  i.e.  parallel  2D  flow.  The  perturbation  stream  function  representing 
a  single  Tollmien-Schlichting  wave  is  then  assumed  to  be  of  the  form: 

'J‘(x,z,t)  =  Q(z)eia(-x-m)  (1) 

The  components  of  the  perturbation  velocity  can  then  be  expressed  as: 

w'=3-  =  (p'izW^x-^+cc,  (2) 

oz 

3\t/ 

vv'  =  —  =  -iaQizy^-^+cc.  (3) 

OX 

Herein,  a  and  ft)  are  in  general  complex, 

a  =  ar  +  ioci ,  (4) 

co  =  cor  +  ia>i ,  (5) 

with  ar  being  the  wave  number  (wavelenght  A  =  In/ ar ),  cor  the  circular  frequency  and  a/,  co*  amplification 
rates.  It  can  then  be  distinguished  between  a  spatial  and  a  temporal  theory:  In  the  spatial  theorie  co  is  real, 
and  amplification  rates  Of/  <  0  denote  spatial  amplified  disturbances  in  flow  direction,  while  in  the  temporal 
theorie  a  is  real  and  amplifications  in  time  are  given  for  amplification  rates  ft)/  >  0  and  given  positions  x. 

Using  the  spatial  theory  and  introducing  the  perturbation  velocities,  eq.  2  and  3,  into  the  Navier-Stokes 
equations,  one  can  obtain  the  Orr-Sommerfeld  equation  for  the  amplitude  <p: 

cp""  -  2 a2cp"  +  a4(p  =  iRe[(aU  -  co)(cp"  -  a2<p)  -  aU"cp\  .  (6) 

Using  homogenous  boundary  conditions  for  the  perturbations  at  the  wall  and  for  large  values  of  y,  the 
problem  of  stability  is  now  reduced  to  an  eigenvalue  problem  of  eq.  6,  assuming  that  the  mean  flow,  Reynolds 
number  and  circular  frequency  ft)  are  given.  Then,  the  solution  yields  the  most  unstable  complex  eigenvalue 
a  =  ar  +  ioci. 

2.1.2  Transition  Prediction  for  Unsteady  Flow  with  the  e^-Method 

A  quite  successful  method  for  transition  predictions  was  derived  from  observations,  that  the  location  of  the 
final  transition  phase  (break  down)  is  in  many  cases  dominated  by  the  behaviour  of  the  primary  instabilities 
with  exponential  growth.  It  has  been  found  that  the  point  where  the  boundary  layer  becomes  fully  turbulent 
correlates  strongly  with  a  certain  amplification  factor  of  the  most  unstable  primary  wave  that  is  calculated 
from  the  point  of  neutral  stability  up  to  the  location  with  fully  turbulent  flow.  These  findings  constitute  the 
so  called  ^-method. 

The  assumptions  needed  to  apply  the  method  to  the  problem  of  predicting  transition  of  LSBs  include: 

•  Initial  and  external  disturbances  of  the  laminar  boundary  layer  are  small.  These  are  the  surface  rough¬ 
ness,  external  turbulence,  acoustic  disturbances,  noise,  and  probably  others. 
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ORGANIZATION 


•  The  location  of  the  final  transition  phase  (break  down)  is  governed  by  the  behaviour  of  primary  20- 
waves  (Tollmien-Schlichting  waves). 

•  The  laminar  boundary  layer  is  thin  and  grows  slowly  with  x. 

The  last  assumptions  may  appear  questionable,  at  first  sight.  However,  extensive  numerical  experience  from 
numerical  stability  computations  and  comparisons  with  the  results  of  DNS  by  U.  Rist  and  his  co-workers 
indicates  that  the  assumption  of  parallel  flow  needed  for  the  LST  formulated  above  is  valid.  That  is,  the 
growth  rates  of  primary  instabilities  computed  by  LST  compares  very  well  with  the  rates  predicted  by  well 
resolved  DNS  [10],  [11]. 

In  the  spatial  theory,  the  overall  amplification  factor  A(x)/Ao  of  the  perturbation  amplitude  A  is  built  by 
an  integration  of  the  local  amplification  rates  af. 


A(x)  ~J<*idx 

v  J  e  xo 

^0 


(7) 


This  integration  is  done  for  a  fixed  number  of  values  for  the  circular  frequency  co.  The  so  called  N-factor  is 
then  obtained  by  taking  the  maximum  value  of  the  amplitude  exponent 


N  —  max 


(8) 


and  thus  forms  an  envelope  over  the  investigated  modes.  This  is  shown  in  Fig.  2  for  the  example  of  of  the 
flow  over  a  flat  plate. 


Fig.  2:  Application  for  flat  plate,  from  D.Arnal  [1  ] 


The  energy  of  a  wave  propagates  with  the  group  velocity 

dcor 
V8  =  Ja’ 


(9) 


which  gives  the  link  between  the  spatial  and  temporal  amplification  rates,  known  as  Gaster  transformation: 


Cdj  —  a/v, 


(10) 
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Then,  while  the  energy  of  a  mode  is  convected  downstream  to  a  position  x  =  x  +  vgA t,  its  perturbation 
amplitude  is  amplified  by 


A  n  =  In 


M*) 

A(x) 


t-\-At 

J  (Oi(x)dt 

t 


(11) 


This  can  be  used  to  calculate  amplitude  exponents  for  unsteady  flows  using  a  simple  temporal  integration 
scheme.  For  the  new  position  x  and  for  a  single  mode  at  the  time  t  +  At  this  yields: 


n(x)  —  n(x)  +  A n(x)  . 


(12) 


The  N-Factor  is  than  again  obtained  by  taking  the  maximum  value  of  the  amplitude  exponents  n  of  the 
single  modes. 

The  transition  location  is  found  by  assuming  a  critical  N-factor  at  which  transition  is  completed  and 
comparing  it  to  the  calculated  N-Factor  distribution.  For  external  flows  over  airfoils  and  wings  in  wind 
tunnels  the  main  factor  influencing  the  critical  N-factor  is  the  free-stream  turbulence  level  of  the  test  section. 
An  empirical  correlation  between  the  turbulence  level  Tu  and  the  critical  N-factor  Ncr[t  is  given  by  Mack 
[13]: 


Ncrit  =  -8A3-2Aln(Tu) 


(13) 


2.1.3  Approximate  Envelope  Method  for  Steady  Flow 

Exact  solutions  of  the  two-dimensional  boundary  layer  equations  for  laminar  flow  can  be  obtained  by  solving 
the  Falkner-Skan  equation 


/"'  +  /"  +  &(l-/2)  =  0,  (14) 

where  the  non-dimensional  streamfunction  /(rj)  depends  only  on  the  wall-normal  coordinate.  The  solution 
describes  the  flow  past  a  wedge  with  the  wedge  angle  equal  to  7 r/3/*  and  consist  of  self-similar  velocity 
profiles.  These  were  studied  in  detail  by  D.R.  Hartree,  therefore  the  parameter  fa  is  referred  to  as  Hartree- 
parameter  and  the  corresponding  velocity  profiles  as  Hartree  profiles.  For  every  Hartree  profile  a  distinct 
shape  parameter  H  can  be  affiliated. 

Each  Hartree  profile  can  now  be  analysed  similar  to  the  flat  plate  as  shown  in  Fig.  2,  which  is  the  special 
case  of  the  flow  past  a  wedge  with  an  angle  of  0°  (j fa  =  0).  This  yields  envelope  curves  similar  to  the  one 
in  Fig.  2,  in  dependence  of  the  shape  parameter  H.  These  curves  can  be  approximated  by  linear  functions. 
These  functions  are  governed  by  two  parameters,  one  for  the  intersection  with  the  abscissa  and  one  for  the 
gradient.  Polynomials  may  be  defined  that  yield  numerical  values  of  these  parameters  and  these  polynomials 
are  meant  when  referring  to  “the”  envelope.  The  envelope  used  for  the  calculations  in  the  present  paper  has 
been  created  by  Wiirz  [6]  using  tables  calculated  by  Arnal.  In  order  to  apply  this  envelope  for  the  flow  around 
an  airfoil  it  is  assumed  that  the  velocity  profiles  of  the  airfoil  boundary  layer  are  similar  to  the  Falkner-Skan 
flows  with  the  same  shape  parameter  H.  Therefore,  they  are  assumed  to  exhibit  the  same  stability  behaviour. 


2.2  Turbulence  Modelling 

Turbulence  modelling  for  engineering  applications  is  usually  based  on  the  Reynolds  averaged  Navier-Stokes 
equations.  In  this  approach  the  contribution  of  turbulence  to  the  mean  momentum  is  accounted  for  by  the 
divergence  of  the  so-called  Reynolds  stress  tensor,  which  has  to  be  modelled  in  terms  of  mean  flow  quantities. 

Most  models  rely  on  the  Boussinesq  hypothesis,  namely  that  the  Reynolds  stress  tensor  is  proportional  to 
the  mean  (traceless)  shear  tensor,  i.  e. 


pRij  =  -iit 


dUi  dUj 
dxj  dxi 


2dUk 
3  dxk 


2 

+  -pkSij, 


(15) 
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where  U[  are  the  mean  velocity  components  and  k  is  the  kinetic  turbulence  energy.  The  proportionality 
coefficient  pt  is  the  eddy  viscosity  and  must  be  provided  by  the  respective  model. 

In  algebraic  turbulence  models,  like  that  of  Baldwin  and  Lomax  [14],  the  eddy  viscosity  is  provided  by 
some  algebraic  relation,  often  based  on  Prandtl’s  mixing  length  concept. 

Improvements  are  achieved  using  one  or  two-equation  models.  The  one-equation  model  of  Spalart-Allmaras 
[15]  uses  a  transport  equation  for  a  value  strongly  related  to  the  turbulent  viscosity  with  rather  heuristic  mo¬ 
delling  of  turbulence  diffusion,  production  and  destruction.  However,  this  model  is  numerically  very  robust 
and  therefore  widely  used  for  applications.  Two-equation  models  exploit  transport  equations  for  the  kinetic 
turbulence  energy  k  and  some  length  scale  of  the  turbulent  eddies,  e.  g.  £  or  co.  A  model  of  this  class  that  is 
becoming  increasingly  popular  is  the  SST  model  of  Menter  [16],  which  reads 
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(17) 


This  model  combines  the  advantages  of  k-£  and  k-co  models  by  gradually  changing  the  closure  coefficients 
/3,  /3*,  y,  Gj,  Ok  and  0®  via  a  blending  function  F\ .  Additionally  the  model  has  been  sensitized  to  predict 
separation  more  reliably  by  defining  the  eddy  viscosity  as 


Pt  = 


a\pk 

max^ift);/^)  ’ 


(18) 


where  Cl  is  the  absolute  value  of  the  vorticity  and  F2  is  a  function  depending  on  the  wall  distance.  Neverthe¬ 
less  all  these  models  suffer  from  the  weaknesses  of  the  Boussinesq  assumption  (15). 

This  drawback  is  overcome  by  solving  the  transport  equations  for  the  individual  Reynolds  stress  com¬ 
ponents  R[j,  which  read 


d(PR 


dt 


JJ  ~ckx  ~  pP ;j  +  pDij  —  p£ij  +  pTLij. 


(19) 


In  contrast  to  two-equation  models,  in  this  approach  the  production  term  pPij  =  — pR ik^r  -p%Sisexact’ 
while  the  diffusion  term  pDz/,  the  destruction  term  and,  in  addition,  the  redistribution  term  pYlij  need 

to  be  modelled.  In  this  context  simple  gradient  diffusion  is  assumed,  i.  e.  pDz/  =  ^  ^  /~\  dR> 


C P+CsPe/e )%l 


dxk 

while  for  the  destruction  term  it  is  common  practice  to  assume  isotropy,  i.  e.  £/y  =  \z8ij.  Like  in  two-equation 
models,  £  (or  co)  must  be  provided  by  an  additional  length  scale  equation. 

A  redisribution  model  that  performs  very  well  in  homogenous  turbulence  is  the  SSG  model  [17],  reading 

n  ij  =  —  (cj  ^ £  +  ~c[  ^Pkk'j  bij  +  C^£  i^bikkkj  ~  ^ bmnbmn8i +  ^3  —  C3  \J  bmnbm}^j  kSij 

+C4k  (bikSjk  +  bjksik  -  -bmnSmn8ij  \  -  C5k  (bikWkj  -  Wikbkj)  (20) 


In  this  equation  bij  =  Rij/2k  —  5/y/3  is  called  the  anisotropy  tensor,  and  S[j  and  W[j  represent  the  shear  and 
the  rotation  tensors,  respectively.  For  the  wall  bounded  flows  investigated  here  the  model  has  been  combined 
with  the  somewhat  simpler  stress-ft)  model  of  Wilcox  [18]  near  walls  by  gradually  changing  the  coefficients 
C{  via  the  blending  function  F\  of  the  SST  model  [19].  Equation  (17)  is  used  for  providing  the  length  scale. 
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3  NUMERICAL  METHODS 

3.1  Navier  Stokes  Code  FLOWer 

The  Navier-Stokes  solver  FLOWer  [2<  ]  uses  a  block-structured  computational  domain  around  the  aerody¬ 
namic  configuration.  The  flow  equations  are  discretised  based  on  the  finite-volume  approach,  either  with 
the  cell-vertex  or  the  cell-centred  formulations.  Various  options  are  implemented  for  treating  the  convective 
fluxes.  These  are  central  schemes  with  scalar  or  matrix  valued  dissipation  and  upwind  schemes  as  well.  Inte¬ 
gration  in  time  is  performed  by  explicit  multistage  schemes.  Implicit  residual  averaging,  multigrid  and  local 
time  steps  improve  convergence  for  steady- state  flow  computation.  Time-accurate  computations  are  suppor¬ 
ted  by  using  the  dual  time  stepping  approach.  Both  convergence  and  accuracy  of  low-speed  flow  simulations 
are  enhanced  by  preconditioning.  The  code  is  optimised  for  vector  computers.  Parallel  computations  are  ba¬ 
sed  on  MPI  and  the  use  of  a  high-level  communication  library.  A  variety  of  turbulence  models  is  implemented 
into  FLOWer.  These  include  algebraic  eddy- viscosity  models,  one-  and  two  equation  models,  algebraic  stress 
models  and  Reynolds  stress  models.  Among  these  the  one-equation  model  of  Spalart-Allmaras  and  the  k—  CO 
model  of  Wilcox  are  standard  options,  but  their  variants  SST,  LEA,  EARSM  have  shown  improved  results 
for  specific  applications,  particularly  with  turbulent  flow  separations. 

3.2  Linear  Stability  Equations  Solver  Coast3 

The  numerical  method  used  for  determining  local  amplification  of  disturbances  is  more  general  compared  to 
the  stability  theory  based  on  the  Orr-Sommerfeld  equation  as  presented  in  chapter  2.1.  Coast3  treats  laminar, 
compressible  boundary  layers  [21  ].  The  boundary  layer  is  assumed  to  be  a  parallel  flow.  The  harmonic  wave 
ansatz  is  applied  to  the  variables  u,v,w,p,T.  It  results  in  a  system  of  five  linear  differential  equations  with  an 
option  to  include  curvature  effects.  Coast3  solves  the  eigenvalue  problem  for  the  complex  eigenvalue  0),  whe¬ 
re  the  wave  numbers  a  and  /3  (for  3D-Cases)  are  prescribed  by  the  user.  Hence  the  temporal  stability  problem 
is  solved.  The  system  of  differential  equations  is  discretised  with  symmetrical  second-order  differences.  This 
yields  a  5n-dimensional  complex  band  matrix  where  n  is  the  number  of  grid  points  in  wall-normal  direction. 
The  numerical  eigenvalue  computation  is  accomplished  with  a  generalised  inverse  Rayleigh  iteration  [22] 
that  takes  into  account  the  banded  structure  of  the  algebraic  problem.  The  code  provides  efficient  searching 
strategies  that  can  be  used  to  find  the  range  of  unstable  modes  for  a  given  flow  problem.  A  data  base  method 
to  find  the  range  of  amplified  spanwise  wave  numbers  is  offered  as  well.  Also  an  option  for  boundary  layer 
computations  for  a  given  pressure  distribution  is  available. 

4  IMPLEMENTATION  ISSUES 

4.1  Extraction  of  Boundary  Layer  Parameters 

In  order  to  calculate  N-factor  distributions  with  the  approximate  envelope  method  as  well  as  the  linear  stabi¬ 
lity  solver,  several  integral  boundary  layer  parameters  have  to  be  extracted  from  the  Navier-Stokes  solution. 
Therefore,  the  boundary  layer  edge  has  to  be  detected  first.  As  the  pressure  can  be  assumed  to  be  constant 
normal  to  the  boundary  layer,  the  boundary  layer  edge  velocity  Ue  can  be  calculated  from  the  wall  pressure, 
pw ,  using  the  compressible  Bernoulli  equation: 


Ue  = 


u 2 

^  oo 


2  K  poo 

K-  1  Poo 


(21) 


Then,  the  boundary  layer  thickness  8  is  found  at  the  point  in  wall-normal  direction  where  a  specified  fraction 
U  =  c  •  Ue  of  the  edge  velocity  is  reached,  with  c  —  0.99  as  a  typical  value.  If  U  is  not  reached,  then  the 
point  with  the  maximum  velocity  is  taken  instead.  These  criterions  do  not  perform  well  close  to  the  leading 
stagnation  point.  Therefore,  all  points  with  a  surface  pressure  coefficient  cp  of  more  than  95  percent  of 
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the  stagnation  value  are  excluded.  This  can  be  safely  done  because  disturbances  are  damped  in  this  region 
anyway. 

As  the  wall  boundary  layer  has,  in  contrast  to  the  parallel  flow  assumption  of  the  boundary  layer  theory,  a 
small  velocity  component  v  normal  to  the  wall,  the  absolute  velocity  is  used  for  Ue  here, 

Ue  =  \/m2  +  v2,  (22) 

while  the  wall  tangential  velocity  component  is  used  for  U  to  calculate  the  boundary  layer  parameters: 


ORGANIZATION 


U  =  utan. 


(23) 


The  displacement  thickness  <5*  and  the  momentum  loss  thickness  6  can  be  determined  by  integration  from 
the  wall  to  the  boundary  layer  edge: 


<5* 


(24) 


Using  the  above  values  for  5*  and  0,  the  shape  parameter  and  the  Reynolds  number  based  on  6,  the  velocity 
Ue  and  the  viscosity  ve  at  the  boundary  layer  edge  can  be  determined: 


<5* 

H=  y  ,  Ree 


Ue6 


Ve 


(25) 


4.2  Treatment  of  Stability  Data  Close  to  the  Numerical  Transition  Location 

The  numerical  update  of  the  transition  location  requires  special  attention:  Switching  on  the  source  terms  of 
the  turbulence  model  at  the  predicted  transition  location  will  numerically  affect  the  laminar,  boundary  layer 
just  upstream  of  that  point,  by  the  inherent  discretisation  errors  of  the  numerical  scheme.  This  numerical  error 
can  have  adverse  effects  on  the  further  transition  prediction  result  if  no  adequate  treatment  is  introduced  into 
the  transition  location  update.  In  our  experience,  for  the  LST-solver  as  well  as  the  approximative  envelope 
methods,  this  is  mainly  an  issue  in  attached  flows.  Also,  separated  flows  are  effected  if  turbulence  models  are 
used  which  react  abruptly  to  their  activation,  like  algebraic  turbulence  models.  We  have  therefore  developed 
a  new  method  for  extracting  reliable  transition  location  estimates  from  stability  results  based  on  RANS 
mean  flow  data.  The  method  is  based  on  stream  wise  linear  extrapolation  of  N-factors  from  flow  regions  with 
reliable  laminar  mean  flow  into  the  transitional  or  even  turbulent  locations.  The  method  searches  for  non¬ 
physical  gradients  of  the  N-factor  or  local  maxima  that  are  caused  by  the  numerical  discretisation  errors  close 
to  transition  in  order  to  find  the  extrapolation  region  in  an  automated  way.  A  sample  using  the  approximative 
envelope  method  is  shown  in  Fig.  3.  The  result  is  taken  from  an  airfoil  computation  where  transition  takes 
place  above  a  LSB.  Here,  the  investigation  of  the  turbulent  region  behind  as  well  as  the  influenced  region 
ahead  of  the  transition  location  results  in  a  strong,  nonphysical  increase  of  the  N-Factor,  which  is  overcome 
by  the  linear  extrapolation.  Without  the  extrapolation,  not  only  the  transition  location  would  be  false,  but 
also  oscillations  could  occur,  as  sometimes  the  local  maximum  ahead  of  the  strong  increase  of  the  N-Factor 
is  detected  at  transition  location.  Much  stronger  oscillations  can  appear  using  the  linear  stability  solver  in 
attached  flows,  as  shown  in  Fig.  4.  Here,  with  increasing  time  it  can  happen  that  no  transition  location  is 
found  at  all.  As  the  transition  location  is  then  set  to  the  trailing  edge,  one  obtains  a  strongly  amplified  laminar 
region  and  thus  the  transition  location  moves  to  the  previous  position.  The  linear  extrapolation  method  can 
therefore  be  used  to  avoid  these  unphysical  oscillations.  Note  that  the  e^-method  itself  is  an  extrapolation 
technique  of  linear  stability  results  to  the  transition  point  upstream  of  which  nonlinear  effects  dominate  in 
reality.  Therefore,  a  local  extrapolation  of  N-factors  is  an  acceptable  modification  of  the  overall  e^-method. 
Finally,  under-relaxation  in  the  iterative  update  of  the  numerical  transition  location  may  be  used  in  order  to 
damp  oscillations  in  the  iterative  transition  location  update. 
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Fig.  3:  N-Factor  distribution  for  separated  flow  Fig.  4:  N-Factor  distribution  for  attached  flow 
using  the  approximative  envelope  method  using  the  linear  stability  solver 


4.3  Transition  Prediction  Procedure  Using  the  Linear  Stability  Solver 

4.3.1  Approach  for  Steady  Flows 

For  the  calculation  of  steady  state  flows  the  linear  stability  solver  is  used  in  usual  ways.  First,  the  frequency 
range  of  amplified  modes  is  estimated.  Therefore  all  eigenvalues  at  an  initial  station  are  searched.  Then,  the 
mode  with  the  most  amplified  eigenvalue  at  this  station  is  tracked  up-  and  downstream  until  it  is  not  amplified 
anymore.  At  these  two  positions,  the  modes  with  the  highest  and  lowest  frequency  that  are  still  amplified  are 
determined.  Second,  a  N-factor  analysis  is  performed  in  the  previously  determined  frequency  range  for  a 
number  of  fixed  frequencies.  For  all  these  given  frequencies  F,  the  most  amplified  modes  are  tracked  up- 
and  downstream  until  they  are  damped.  Their  temporal  amplification  rates  are  calculated  in  this  region  and 
tranformed  in  spatial  ones  using  the  Gaster  transformation.  These  are  then  integrated  to  N-factor  distributions 
in  the  amplified  region  for  every  frequency.  The  envelope  over  all  these  N-factor  distributions  then  yields  the 
N-factor  curve  which  is  used  to  determine  the  transition  location  by  comparsion  with  a  critical  N-factor. 

4.3.2  Approach  for  Unsteady  Flows 

For  a  given  frequency  multiple  modes  with  different  wavenumbers  and  amplification  rates  may  be  found 
by  the  numerical  method.  As  the  N-factor  distribution  for  a  given  frequency  is  now  computed  with  a  time- 
dependant  scheme,  it  has  to  be  made  sure  that  at  every  timestep  the  same  modes  are  analysed.  This  is  not 
guaranteed  with  the  above  approach  for  steady  flows,  which  is  why  it  is  now  used  in  a  modified  way.  As  the 
mean  flow  changes  with  time,  also  the  frequency  range  of  amplified  modes  changes.  Accordingly,  the  inves¬ 
tigated  range  of  frequencies  has  to  be  adapted  during  an  unsteady  calculation.  For  the  present  computations, 
this  is  circumvented  by  selecting  a  fixed  frequency  range  and  checking  afterwards  if  the  frequencies  that 
were  relevant  for  the  transition  were  investigated  at  all  times.  The  N-Factor  analysis  is  then  performed  once 
for  every  mode,  using  the  modes  wavelengths  and  amplification  rates  from  the  last  timestep  as  initial  values 
for  the  solver  to  make  sure  that  the  same  mode  is  investigated  again.  As  new  amplified  modes  can  form 
in  unsteady  flows,  these  are  searched  after  every  timestep  along  the  airfoil  for  the  investigated  frequencies. 
Also,  modes  that  are  damped  to  the  point  where  the  amplitude  exponent  is  zero  over  the  whole  airfoil  are 
identified  and  deleted. 
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5  NUMERICAL  FLOW  COMPUTATIONS 


ORGANIZATION 


Steady  state  flow  solutions  are  often  not  obtained  for  flows  with  LBS's,  even  with  Reynolds  averaging  of  the 
flow  equations.  As  a  consequence,  the  RANS  solver  is  run  in  its  time-accurate  mode  using  the  dual  time¬ 
stepping  approach,  even  for  steady  onset  conditions.  As  a  rule  of  thumb  we  use  a  time  step  length  of  At  « 
0.01,  where  the  time  is  made  non-dimensional  by  using  the  onset  flow  velocity  and  the  chord  length.  Figure 
5  indicates  that  this  time  is  sufficiently  small  to  resolve  the  main  viscous  boundary  layer  interactions.  For 
unsteady  onset  conditions  there  remains  a  very  small  influence  on  the  temporal  distribution  of  the  transition 
location,  as  seen  in  Fig.  6.  As  the  influence  on  the  distribution  of  the  force  coefficients  is  negligible,  no 
shorter  timestep  than  At  «  0.01  was  chosen  to  minimise  the  computational  time. 


Fig.  5:  Variation  of  timestep  length  for  steady  on-  Fig.  6:  Variation  of  timestep  length  for  unsteady 
set  conditions  at  Re  =  1.2  •  106  onset  conditions  at  Re  =  60000 


5.1  Steady  Onset  Conditions 

5.1.1  Flow  Computations  of  XIS40MOD  Airfoil 

Numerical  computations  of  the  research  airfoil  XIS40MOD  are  chosen  because  there  exists  detailed  experi¬ 
mental  data  of  Wiirz  [6]  from  a  low-disturbance  wind  tunnel.  The  measured  data  include  pressure  distribu¬ 
tions  and  boundary-layer  flow  velocities  measured  with  a  hot  wire  (mean  flow  and  flow  disturbances).  The 
airfoil  is  specifically  designed  to  produce  a  large  laminar  separation  bubble  on  its  upper  surface  (see  Fig.  7). 
The  computational  mesh  shown  has  528  cells  in  surface-tangential  direction  whereas  normal  to  the  surface 
96  cells  are  used.  This  is  the  standard  mesh.  The  mesh  points  are  clustered  around  the  leading  edge  and 
around  the  expected  LBS.  About  40  out  of  96  wall-normal  cells  are  located  in  the  laminar  boundary  layer 
(not  shown  here).  More  refined  meshes  have  been  generated  in  order  to  obtain  grid  convergence.  The  fine 
mesh  has  1056x192  cells  and  a  super-fine  mesh  with  2112x384  cells  was  also  used.  The  flow  conditions  are 
chosen  as  Re  =  1 .2  •  106  and  a  =  —3°,  according  to  the  experiment.  The  first  round  of  computations  is  perfor¬ 
med  using  a  fixed  transition  location  at  x/Z  =  0.76  on  the  upper  airfoil  surface  in  order  to  address  the  effect  of 
spatial  0  errors  on  the  numerical  solutions.  The  results  are  displayed  in  Fig.  8  and  Fig.  9.  The  computations 
are  performed  with  the  SST  turbulence  model.  The  pressure  distribution  and  skin  friction  are  well  resolved 
in  the  region  of  attached  flow  using  the  528x96  mesh.  However,  it  takes  a  grid  of  around  1056x192  points 
in  order  to  capture  the  laminar  separation  accurately.  The  grid-converged  computation  exhibits  a  pressure 
rise  in  the  rear  part  of  the  bubble  that  is  not  as  steep  as  in  the  experimental  data.  Also,  the  wall  shear  stress 
downstream  of  reattachment  is  significantly  smaller  in  the  computation  compared  with  the  measurement. 

The  next  computations  are  used  to  investigate  the  effect  of  varying  the  turbulence  model.  These  com¬ 
putations  are  performed  on  the  528x96  grid.  The  results  of  the  algebraic  model  of  Baldwin  and  Lomax  are 
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Fig.  7:  Standard  computational  grid  for  XIS40MOD  airfoil  and  computed  distribution  of  pressure  coefficient  at 

Re  =  1.2-  106,a  =  —3° 


compared  with  the  one-equation  model  of  Spalart  and  Allmaras,  two  variants  of  the  k  —  co  model  and  the  SSG 
Reynolds-stress  model  in  Fig.  10  and  11.  Note  that  different  numerical  transitions  locations  were  used  as  the 
individual  models  predict  a  rather  different  rise  of  turbulence  and  turbulent  eddy  viscosity  downstream  of 
the  nominal  transition  location.  For  example,  it  appears  that  the  growth  of  turbulent  kinetic  energy  predicted 
by  the  original  k  —  CO  model  is  rather  slow,  compared  to  the  others.  Therefore,  in  order  to  compute  similar 
positions  of  the  turbulent  reattachment  for  this  model  one  has  to  specify  the  nominal  transition  position  si¬ 
gnificantly  upstream.  The  results  indicate  that  the  differential-equation  turbulence  models  predict  a  too  slow 
growth  of  turbulence  in  the  separation  bubble.  Accordingly,  the  rise  of  wall  shear  stress  in  the  re-attachment 
region  is  much  too  slow  either.  This  problem  is  more  striking  for  the  original  k—  CO  model  whereas  the 
Spalart- Allmaras,  the  SST  and  the  SSG  models  perform  somewhat  better.  Note  that  the  SST  model  uses  the 
k  —  £  formulation  far  away  from  walls  and  this  seems  to  be  the  reason  for  the  more  rapid  growth  of  turbulence 
compared  to  the  k  —  co  model.  This  holds  true  for  the  SSG  model  as  well.  The  algebraic  Baldwin-Lomax  mo¬ 
del,  on  the  other  hand,  yields  an  instantaneous  switch  to  fully  turbulent  flow  right  at  the  nominal  transition 
point.  It  turns  out  that  this  results  in  rather  close  agreement  of  the  predicted  pressure  distribution.  However, 
all  models  share  the  unsatisfactory  result  that  the  wall  shear  stress  is  significantly  underpredicted  downstream 
of  re-attachment.  Here,  the  high  complexity  of  the  Reynolds-stress  model  seems  have  a  positive  effect,  as 
it  predicts  the  highest  level  of  turbulent  shear  stress  of  all  models.  However,  the  results  of  the  currently  im¬ 
plemented  Reynolds-stress  model  can  only  be  understood  as  step  in  the  right  direction.  Please  note  that  the 
experimental  hot-wire  set  up  could  not  measure  reverse  flow.  Therefore,  the  positive,  experimental  values  of 
Cf  within  the  separated  flow  region  are  false. 

The  flow  structures  of  the  laminar  separation  bubble  predicted  by  the  algebraic  model  and  the  SST  model 
are  displayed  in  Fig.  12.  The  vertical  co-ordinate  is  enlarged  by  a  factor  of  3  in  order  to  make  the  differences 
visible.  The  SST  model  predicts  a  shallow  shape  of  the  separation  bubble  whereas  the  bubble  of  the  algebraic 
model  closes  much  more  rapidly.  A  fast  turning  vortex  is  created  in  the  aft  part  of  the  LSB  by  the  sudden 
onset  of  turbulent  viscosity  with  the  Baldwin-Lomax  model.  This  vortex  corresponds  to  a  local  peak  in  the 
boundary-layer  edge  velocity  at  x/c  =  0.8,  which  is  visible  in  both,  computation  and  experiment. 

Next,  the  effect  of  grid  density  on  the  predicted  N-factor  is  investigated.  Figure  13  displays  the  N-factor 
distribution  that  results  either  from  using  the  approximate  envelope  method  or  from  solving  the  local  linear 
stability  equations  with  the  Coast3  code.  Coast3  is  used  in  two  ways:  In  one  way  only  the  pressure  distri¬ 
bution  from  FLOWer  is  used,  and  the  affiliated  boundary  layer  method  “Coco”  generates  the  velocity  and 
temperature  profiles  up  to  the  point  of  separation.  In  the  other  way,  the  profiles  are  used  directly  from  the 


RTO-MP-AVT-111 


UNCLASSIFIED/UNLIMITED 


3-11 


UNCLASSIFIED/UNLIMITED 


RANS  Simulation  of  the  Transitional  Flow  Around  Airfoils  at 
Low  Reynolds  Numbers  for  Steady  and  Unsteady  Onset  Conditions 


Fig.  8:  Boundary  layer  edge  velocity  for  different  Fig.  9:  Coefficient  of  wall  shear  stress  for  different 
grid  densities  grid  densities 


Fig.  10:  Boundary  layer  edge  velocity  for  different  Fig.  11 :  Coefficient  of  wall  shear  stress  for  diffe- 
turbulence  models  rent  turbulence  models 


Navier-Stokes  method  over  the  whole  airfoil.  As  the  accuracy  of  the  laminar  mean  flow  improves  on  finer 
meshes  one  observes  a  small  rise  of  N-values  with  an  increasing  number  of  points  in  the  boundary  layer.  For 
about  70  points  in  the  boundary  layer  good  agreement  with  the  measurement  as  well  as  the  N-factor  distribu¬ 
tion  based  on  the  pressure  distribution  is  found.  The  approximative  envelope  method  also  shows  quite  good 
agreement  in  the  region  of  attached  flow.  The  increasing  disagreement  above  the  LSB  is  due  to  the  fact  that 
the  method  assumes  Hartree-profiles  (see  2.1.3),  which  is  invalid  in  the  region  of  the  LSB. 

Finally,  the  iterative  convergence  of  transition  location  with  time  is  displayed  in  Fig.  14.  The  computation 
was  started  from  an  initial  transition  location  close  to  the  trailing  edge.  This  provokes  laminar  flow  separation 
in  the  initial  phases  of  the  computation.  The  separation  then  grows  until  the  N-factor  becomes  larger  than 
the  critical  value.  Then,  transition  jumps  upstream  and  re-attachment  is  obtained  by  turbulent  momentum 
transport.  Finally,  the  transition  location  converges  continuously  to  the  steady-state  result.  This  takes  about 
300  physical  time  steps  of  At  =  0.01,  according  to  Fig.  14.  In  conclusion  we  notice  that  the  computations 
of  the  XIS40MOD  airfoil  display  a  strong  sensitivity  of  the  computed  LSB  to  spatial  discretisation  errors. 
However,  extraction  of  sufficiently  accurate  boundary  layer  data  for  transition  prediction  is  feasible.  The 
major  problem  in  physical  modelling  of  this  case  seems  to  be  the  turbulence  model.  Two-equation  models 
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Fig.  12:  Comparison  of  flow  streamlines  for  two 
turbulence  models 


Fig.  13:  Sensitivity  of  the  N-factor  distribution  on 
methods  and  grid  density,  SST-model  with 

xtran  =  0.76 


do  not  predict  the  strong  rise  of  shear  stress  that  would  bring  the  pressure  distribution  and  wall  shear  stress 
close  to  the  experimental  observations.  The  good  agreement  of  the  algebraic  model  in  terms  of  pressure 
distribution  seems  fortuitous  as  this  model  simply  neglects  any  history  effects  of  turbulence  associated  with 
convection  and  diffusion. 


5.1.2  Flow  Computations  of  SD7003  Airfoil 

The  SD7003  airfoil  is  designed  for  model  airplanes  that  operate  in  the  Reynolds  number  range  of  0.5  •  105  — 
2  •  105.  The  pressure  distribution  exhibits  a  carefully  controlled  adverse  pressure  gradient  ahead  of  laminar 
separation.  This  strategy  yields  growing  disturbances  ahead  of  the  separation.  Furthermore,  the  adverse  pres¬ 
sure  gradient  in  the  rear  of  the  airfoil  is  rather  limited,  due  to  the  small  airfoil  thickness  of  around  8  percent 
and  its  forward  location  of  the  point  of  maximum  thickness.  The  resulting  separation  bubble  is  therefore  thin 
over  the  full  angle  of  attack  range,  even  for  Reynolds  numbers  below  105.  Hence,  the  airfoil  is  not  suspect 
to  the  break  up  of  the  laminar  separation  bubble  within  its  design  range  and  it  is  a  good  test  case  to  assess 
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computational  capabilities. 

Recently,  2D  phase  locked  PIV  measurements  were  conducted  on  the  upper  side  of  the  SD7003  in  a 
low  turbulence  ( Tu  —  0.0012)  windtunnel  (LNB)  at  Braunschweig  Technical  University.  The  wing  section 
with  chord  length  200 mm  and  span  400mm  was  made  out  of  glass  fibre  reinforced  plastics  and  fitted  in  the 
400mm  x  600mm  x  1200mm  test  section  of  the  windtunnel.  The  measurements  were  carried  out  at  an  angle 
of  attack  of  a  =  4°  and  a  Reynolds  number  of  Re  —  60000  with  the  commercial  PIV  system  Flow  Master  3S 
and  a  frequency  doubled  Quantel  Twin  double  puls  Nd:  YAG-laser.  To  achieve  highly  detailed  measurements 
along  the  upper  side,  the  observation  domain  was  translated  to  9  different  positions,  each  with  a  field  of  view 
of  27.27mm  x  21.81mm.  At  each  of  the  9  locations  1000  statistically  independent  PIV  images  were  recorded 
in  order  to  calculate  the  mean  velocity  and  three  components  of  the  Reynolds  stress  tensor.  For  the  evaluation 
of  all  images  the  commercial  PIV  software  DaVis  6.2  was  used.  Because  of  a  high  seeding  density  it  was 
possible  to  choose  small  interrogation  windows  with  a  size  of  16 pxx  16 px.  Together  with  an  overlap  of  50 
percent  a  spatial  resolution  of  0.17mm  could  be  achieved.  The  absolute  particle  image  displacement  range 
was  determined  to  be  between  0.004 px  and  14.36 px  which  is  in  an  appropriate  range  for  small  measurement 
errors. 

The  computational  standard  mesh  for  the  numerical  investigation  of  the  SD7003  airfoil  contains  528x96 
cells,  as  for  the  previous  case.  Grid  convergence  studies  using  a  coarse  264x48  mesh  in  addition  indicate  that 
this  case  is,  because  of  the  now  lower  Reynolds  number,  less  sensitive  to  grid  density  as  the  XIS40MOD 
airfoil.  This  is  displayed  for  the  flow  at  a  moderate  angle  of  attack  by  using  the  SST  turbulence  model  in 
Fig.  15.  The  distribution  of  wall-shear-stress  coefficient  indicates  that  a  mesh  with  528x96  cells  is  sufficient 
for  LSB  computations  at  low  Reynolds  numbers.  The  non-dimensional  time  step  was  chosen  as  At  =  0.01 
and  around  1000  time  steps  were  needed  for  convergence.  Note  that  the  number  of  inner  iterations  in  the 
dual  time  stepping  method  was  around  20.  Typical  pressure  distributions  for  a  variety  of  angles  of  attack  are 
displayed  in  Fig.  16.  The  LSBs  are  seen  as  local  plateaus.  The  lower  surface  is  fully  laminar  for  all  angles 
larger  than  zero.  The  LSB  moves  from  the  trailing  edge  to  the  leading  edge  as  the  angle  of  attack  is  increased. 


Fig.  15:  Grid  convergence  of  wall  shear  stress  for 
a  =  6°,  SST-model 


Fig.  16:  Computed  pressure  distributions  with  the 
Baldwin-Lomax  model 


Figure  17  compares  the  turbulent  shear  stress  from  the  PIV  measurement  at  a  =  4°  with  the  computa¬ 
tions  for  different  turbulent  models,  using  the  coupled  linear  stability  solver  coast3.  The  transition  location 
in  the  measurement  is  found  extrapolating  the  wedge  that  is  formed  by  the  turbulent  shear  stress  into  the 
upstream  direction.  The  so  found  transition  location  lies  at  about  the  thickest  point  of  the  LSB,  as  the  ef¬ 
fect  of  the  turbulence  is  the  closure  of  the  bubble,  and  can  be  given  to  about  xtran.meas.  =  0.56.  From  that 
point  on  the  turbulent  shear  stress  continuously  grows  up  to  a  maximum  dimensionless  value  of  -0.026, 
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which  is  reached  about  9.5  percent  of  the  chord  length  c  downstream.  Using  a  critical  N-factor  of  Ncr[t  =  8 
which  fits  the  turbulence  level  of  the  windtunnel  according  to  Mack's  relation  from  chapter  2.1.2  and  the 
algebraic  Baldwin-Lomax  (BL)  turbulence  model,  the  transition  location  is  computed  at  xtran,BL  =  0.53.  As 
the  Baldwin-Lomax  model  produces  large  amounts  of  turbulence  instantaneously  at  the  transition  point,  the 
bubble  is  closed  at  about  the  same  point  as  in  the  measurement,  even  though  only  a  nondimensional  shear 
stress  of  about  half  the  measured  value  is  reached.  The  fact  that  the  transition  location  was  calculated  a  few 
percent  upstream  of  the  measured  indicates  that  the  N-Factor  was  chosen  slightly  to  small,  but  this  results  in 
a  good  overall  simulation  with  the  BL-model.  With  the  SST-model,  a  transition  location  of  xtran.ssT  =  0.42 
is  found,  while  the  first  maximum  of  the  nondimensional  shear  stress  (for  the  shown  instantaneous  flow 
state)  is  reached  10.5  percent  of  c  downstream  of  that.  Although  the  unsteady  flow  behavior  given  by  the 
SST-model  appears  much  closer  to  the  experimental  findings  compared  to  the  BL-model,  the  produced  tur¬ 
bulent  shear  stress  is  (at  the  first  maximum)  again  only  half  of  the  measured  value.  As  this  is  not  sufficient 
to  close  the  bubble,  it  becomes  larger  compared  to  the  BL-model  and  sheds  vortices  at  its  end.  Because  the 
thicker  upstream  region  of  the  bubble  yields  larger  amplifications  of  disturbances,  the  transition  location  is 
computed  11  percent  upstream  compared  to  the  BL-model.  The  SSG-model  experiences  the  same  problems 
as  the  SST-model,  but  even  more  pronounced.  Note  that  this  model  shows  the  best  overall  performance  at 
higher  Reynolds  numbers  [19]  and  was  still  the  superior  transport  equation  turbulence  model  at  the  previous 
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Fig.  17:  Measured  and  computed  turbulent  shear  stress  at  a  =  4° 
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case  at  an  already  low  Reynolds  number  of  Re  —  1.2E6.  It  can  be  concluded  that  the  algebraic  BL-model 
is  best  suited  for  the  simulation  of  this  case,  despite  its  unphysical  behaviour.  Large  improvements  could  be 
expected  from  a  model  based  on  transport  equations  which  would  be  able  to  predict  the  strong  rise  of  the 
turbulent  shear  stress  behind  the  LSB. 

5.2  Unsteady  Onset  Conditions 

5.2.1  Flow  Computations  of  SD7003  Airfoil 

To  validate  the  later  described  numerical  results  measurements  were  conducted  in  the  watertunnel  at  Braun¬ 
schweig  Technical  University.  The  facility  has  a  free  stream  turbulence  level  of  Tu  =  0.01,  which  translates 
into  a  critical  N-factor  of  Ncrit  =  2  using  eq.  13.  The  SD7003  wing  section  with  200 mm  chord  length  and 
a  span  of  248mm  was  mounted  on  a  plunging  mechanism,  as  described  in  [23].  The  mean  angle  of  attack 
was  set  to  ao  =  5.5°.  A  Reynolds  number  of  Re  —  60000  resulted  from  the  onset  velocity  at  infinity  of 
Uoo  =  0.3m /s  and  a  water  temperature  of  20°.  As  a  flapping  frequency  of  /  =  0.25 Hz  and  a  plunging  ampli¬ 
tude  of  zi  =  10mm  were  used,  the  resulting  quasisteady  equivalent  angle  of  attack  amplitude  from  plunging 
can  be  calculated  to  £ \  =  =  3°.  The  resulting  effective  angle  of  attack  for  the  airfoil  motion  then  reads 

ae ff  =  OCo  +  C,\sin(2nft)  —  5.5°  +  3°sin(2nft).  The  airfoil  was  driven  by  an  eccentric  mechanism  with  1  per¬ 
cent  derivation  to  the  ideal  sinusoidal  stroke.  Almost  the  same  2D  PIV  set  up  was  used  as  for  the  windtunnel 
measurements  described  above.  The  field  of  view  was  only  slightly  bigger  (31.24mm  x  24.99mm),  which  led 
to  a  spatial  resolution  of  0.2mm.  To  cover  most  of  the  separation  bubble,  two  slightly  overlapping  locations 
of  30mm  distance  were  chosen.  Phase  locked  measurements  were  performed  for  the  four  motion  positions 
top  dead  centre  (TDC)  (aeff  =  5.5°),  downstroke  (aeff  =  8.5°),  bottom  dead  centre  (BDC)  (aeff  =  5.5°)  and 
upstroke  (aeff  =  2.5°).  Additional  measurements  were  carried  out  for  these  four  angles  of  attack  for  steady 
onset  conditions. 

The  following  numerical  results  were  obtained  using  the  e^-method  with  the  extension  for  unsteady  flows 
as  described  in  chapter  2.1.  A  first  verification  of  the  method  was  achieved  by  a  calculation  of  a  case  with 
steady  onset  condition,  first  with  the  steady  N-Factor  integration  schema  and  second,  as  a  restart,  with  the 
unsteady  one.  There,  the  obtained  N-Factor  distribution  converged  to  a  steady  state  solution  that  agreed 
within  plotting  accuracy. 

As  the  BL  turbulence-model  showed  the  best  performance  for  the  steady  onset  condition  calculations  in 
the  previous  case,  it  was  chosen  for  the  first  unsteady  calculation,  which  is  presented  in  the  following.  To  take 
into  account  the  wall  effect  from  the  above  watertunnel  experiment,  a  simple  wall  correction  for  steady  flows 
from  [24]  was  applied.  This  results  into  an  effective  angle  of  attack  of  aeff  =  6°  +  3.3°sin(2nft )  which  was 
used  for  the  calculations.  As  the  calculation  has  to  be  done  for  nondimensional  values,  a  reduced  frequency 
of  k  =  nfc/Uoo  —  0.52  and  a  plunge  amplitude  of  zi/c  —  0.0554  were  used.  The  used  computational  mesh 
is  similar  to  the  previous  case  and  has  a  density  of  472x96  points. 

Figure  18  shows  the  obtained  temporal  distribution  of  transition  location  and  lift  coefficient,  using  a 
coarser  mesh  of  only  236x48  points  for  the  first  three  motion  periods.  The  nondimensional  timestep  was 
chosen  ten  times  too  large  with  At  =  0.085  during  the  first  five  motion  periods,  which  led  to  some  small 
oscillations  in  the  transition  location  during  the  fourth  and  fifth  motion  period.  The  sixth  and  seventh  motion 
periods  were  then  calculated  using  2000  timesteps  with  At  =  0.0085  each.  As  no  significant  differences  are 
found  between  the  end  of  the  sixth  and  seventh  period,  a  periodic  solution  is  found  for  the  seventh  motion 
period. 

Figure  20  shows  the  turbulent  shear  stress  u'v'  obtained  by  a  phase  locked  measurement  during  the  down- 
stroke  of  the  airfoil.  Even  after  averaging  over  1000  frames  there  are  regions  found  with  elevated  shear  stress, 
which  represent  vortices  shed  at  the  end  of  the  LSB.  This  indicates  that  the  periodic  shedding  of  vortices  be¬ 
hind  the  LSB  is  coupled  to  the  plunge  motion,  which  has  a  much  lower  frequency.  The  transition  locations 
for  these  phase  locked  measurements  have  been  determined  similar  to  the  measurement  at  steady  onset  con¬ 
ditions  in  the  windtunnel  described  earlier.  Here,  as  a  LSB  with  shed  vortices  has  been  measured,  the  first 
local  maximum  of  the  bubble  corresponds  to  the  transition  location  according  to  the  onset  of  the  turbulent 
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Fig.  18:  Temporal  distribution  of  transition  locati-  Fig.  19:  Comparison  of  measured  and  calculated 
on  and  lift  coefficient  transition  locations 


shear  stress. 


Fig.  20:  Measured  turbulent  shear  stress  during  downstroke 


The  transition  locations  found  in  the  computations,  in  the  phase  locked  measurements  and  in  the  steady 
state  flow  measurements  with  the  quasisteady  equivalent  angles  of  attack  are  shown  in  Fig.  19.  The  com¬ 
parison  between  the  measurements  and  computations  for  steady  onset  conditions  show  a  large  discrepancy 
between  the  transition  locations.  It  can  be  concluded  that  either  the  water  tunnels  turbulence  level  is  less  than 
measured,  or  the  correlation  of  Mack  is  not  working  properly  for  the  quite  high  turbulence  level  of  Tu  —  0.01 
of  the  watertunnel.  For  the  unsteady  computation,  the  too  low  chosen  N-factor  of  N  =  2  again  results  in  a 
transition  location  upstream  of  the  measurement,  but  now  with  a  smaller  and  about  constant  distance.  It  is 
interesting  to  see  here  that  the  unsteady  transition  locations  (calculation  and  measurement)  are  downstream 
of  the  steady  ones  at  £  =  8.5°,  while  it  is  the  other  way  around  at  the  lower  angles.  This  behaviour  still  needs 
to  be  explained.  The  computed  unsteady  transition  location  shows  a  phase  shift  of  25°  during  the  downstroke 
and  40°  during  upstroke,  which  is  about  the  expected  effect  at  the  examined  reduced  frequency  of  k  =  0.52. 
The  lift  coefficient  however  seems  to  be  more  dependant  on  the  mean  flow  and  does  not  show  much  phase 
shift,  as  seen  in  Fig.  18. 

The  measured  and  computed  LSBs  are  finally  compared  at  the  dead  centres,  downstroke  and  upstroke 
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in  Fig.  21,  together  with  the  nondimensional  total  velocity  U  =  Vu2  +  V2.  In  order  to  be  able  to  draw 
streamlines,  the  v-component  of  the  velocity  was  corrected  with  the  plunging  velocity  of  the  airfoil,  leading 
to  a  velocity  of  zero  on  the  surface  of  the  airfoil.  It  can  be  seen  that  despite  the  same  quasisteady  angle  of 
attack  of  C  =  5.5°,  there  are  completely  different  flow  states  found  at  TDC  and  BDC.  While  the  LSB  that 
formed  during  the  downstroke  is  clearly  visible  at  BDC,  it  is  vanishing  during  upstroke  and  can  hardly  be  seen 
at  all  at  TDC.  A  transition  location  can  not  be  determined  in  the  measurement  for  TDC.  During  downstroke, 
the  calculated  LSB  is  clearly  higher  than  the  measured  one.  Correspondingly,  the  point  of  laminar  separation 
Xsep/c  =  0.05  is  situated  a  little  ahead  of  the  measured  one  at  xsep/c  =  0.065.  At  BDC  a  vortex  is  shed  from 
the  LSB  in  the  simulation,  which  leads  to  about  the  same  thickness  of  the  LSB  in  the  front  part  as  in  the 
measurement.  The  point  of  laminar  separation  is  again  somewhat  ahead  in  the  calculation. 


ORGANIZATION 


6  SUMMARY 

Transitional  flows  around  airfoils  were  calculated  for  steady  as  well  as  unsteady  onset  conditions  using  a 
Reynolds-averaged  Navier-Stokes  method.  Hereby,  the  transition  location  was  found  using  the  e^-method 
in  combination  with  a  linear  stability  solver.  In  order  to  do  so,  a  new  temporal  integration  schema  for  the 
^-method  was  used.  It  was  pointed  out  that  the  N-Factor  distribution  needs  special  treatment  in  the  vicinity 
of  the  transition  location,  i.e.  extrapolation. 

For  the  flows  with  steady  onset  conditions  it  was  shown  that  high  grid  densities  are  necessary  to  obtain 
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grid  convergence.  It  was  shown  that  the  main  remaining  sensitivity  for  the  flows  with  LSB  is  the  turbulence 
model,  and  that  further  improvement  in  that  field  is  needed. 

A  first  solution  was  presented  for  the  transitional  flow  around  a  moving  airfoil  using  the  new  temporal 
schema.  It  was  found  for  the  investigated  case  that  very  different  flows  form  at  the  quasisteady  identical 
angles  of  attack  at  the  dead  centres,  thus  the  flow  has  an  unsteady  character.  The  computed  unsteady  transition 
locations  were  promising  compared  to  the  measurement. 
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